﻿using System;
using Ewk.Math.ComputationTypes;

namespace Ewk.Math.Algebra.Algorithms.Matrices.Rank
{
    class GaussianElimination<T> : IRankComputer<T>
        where T : struct, IComparable, IFormattable, IConvertible, IComparable<T>, IEquatable<T>
    {
        public int ComputeRank(Matrix<T> matrix)
        {
            if (matrix == null)
            {
                return 0;
            }

            SortRowsToRowEchelonFormBestEffort(matrix);
            CreateRowEchelonForm(matrix);

            var rank = 0;
            for (var i = 0; i < matrix.RowCount; i++)
            {
                var firstColumnNumberContainingNonZero = matrix.GetFirstColumnNumberContainingNonZero(0, i);
                if (firstColumnNumberContainingNonZero == MatrixSearch.NotFound)
                {
                    return rank;
                }

                rank++;
            }

            return rank;
        }

        static void SortRowsToRowEchelonFormBestEffort(Matrix<T> matrix)
        {
            for (var i = 0; i < matrix.ColumnCount; i++)
            {
                for (var j = 0; j < matrix.RowCount; j++)
                {
                    var firstRowNumberContainingZero = matrix.GetFirstRowNumberContainingZero(j, i);
                    if (firstRowNumberContainingZero == MatrixSearch.NotFound) continue;

                    var lastRowNumberContainingNonZero = matrix.GetLastRowNumberContainingNonZero(firstRowNumberContainingZero, i);
                    if (lastRowNumberContainingNonZero == MatrixSearch.NotFound) break;

                    matrix.SwapRows(firstRowNumberContainingZero, lastRowNumberContainingNonZero);
                    j = firstRowNumberContainingZero;
                }
            }
        }

        static void CreateRowEchelonForm(Matrix<T> matrix)
        {
            for (var i = 0; i < matrix.ColumnCount; i++)
            {
                var firstValue = matrix[i, i];
                if (firstValue == Computable<T>.Factory.Zero) continue;

                for (var j = i + 1; j < matrix.RowCount; j++)
                {
                    var secondRowNumberContainingNonZero = matrix.GetFirstRowNumberContainingNonZero(j, i);
                    if (secondRowNumberContainingNonZero == MatrixSearch.NotFound) break;

                    var secondValue = matrix[secondRowNumberContainingNonZero, i];
                    var factor = secondValue/firstValue;

                    for (var k = 0; k < matrix.ColumnCount; k++)
                    {
                        matrix[secondRowNumberContainingNonZero, k] -= matrix[i, k] * factor;
                    }

                    j = secondRowNumberContainingNonZero;
                }
            }
        }
    }
}
